rm(list = ls())
require(stats4)
require(nlWaldTest)
require("multiwayvcov")
require(msm)
require(stargazer)
require(mfx)
require(restriktor)

data21=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december6noon_type7clean.csv", header=TRUE)
data22=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december15_type7clean.csv", header=TRUE)
data23=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_january26_type7(clean).csv", header=TRUE)

dataframe21=data.frame(data21)
dataframe22=data.frame(data22)
dataframe23=data.frame(data23)

#some definitions
expframe=rbind(dataframe21,dataframe22,dataframe23)
expframe=subset(expframe,uid>800)
expframe$correct <-(expframe$state==4&(expframe$response==2|expframe$response==7|expframe$response==8|expframe$response==9))|(expframe$state==3&(expframe$response==1|expframe$response==6|expframe$response==3|expframe$response==5))
expframe$incentive <- (expframe$qid==5)*5+(expframe$qid==1)*40+(expframe$qid==6)*70+(expframe$qid==7)*95
accuracy=c(sum((expframe$incentive==5)*expframe$correct)/sum(expframe$incentive==5),sum((expframe$incentive==40)*expframe$correct)/sum(expframe$incentive==40),sum((expframe$incentive==70)*expframe$correct)/sum(expframe$incentive==70),sum((expframe$incentive==95)*expframe$correct)/sum(expframe$incentive==95))
incentivec=c(5,40,70,95)
uidlist=unique(expframe$uid)

expframe$incentive5 <- (expframe$incentive==5)
expframe$incentive40 <- (expframe$incentive==40)
expframe$incentive70 <- (expframe$incentive==70)
expframe$incentive95 <- (expframe$incentive==95)

expframe$Bresponse <-(expframe$response==2|expframe$response==7|expframe$response==8|expframe$response==9)
expframe$Bstate <- expframe$state==4
expframe$Aresponse <-(expframe$response==1|expframe$response==6|expframe$response==3|expframe$response==5)
expframe$Astate <- expframe$state==3


#check N
#**
length(uidlist)
#**


#Check NIAS and NIAC on the individual level
#NIAS


NIASAframe=data.frame(uid=vector(),A_state_avg_effect=vector(),Pvalue=vector(),incentive=vector())
for (i in 1:length(uidlist)){
for(j in 1:4){
addframe=0
playframe=subset(expframe,uid==uidlist[i] & incentive==incentivec[j])
if(sum(playframe$Aresponse)==0|sum(playframe$Bresponse)==0){
addframe=data.frame(uid=uidlist[i],A_state_coefficient=0,Pvalue=1,incentive=incentivec[j])
NIASAframe=rbind(NIASAframe,addframe)
}else{
model=logitmfx(Aresponse~Astate,data=playframe,robust=TRUE) 
addframe=data.frame(uid=uidlist[i],A_state_coefficient=model[[1]][1],Pvalue=model[[1]][4],incentive=incentivec[j])
NIASAframe=rbind(NIASAframe,addframe)
}
}
}

#Check NIAC 
restrictmat=matrix(c(0,1,0,0, 0,-1,1,0, 0,0,-1,1),nrow=3,ncol=4, byrow=TRUE)

NIACindiframe=data.frame(uid=vector('numeric'),fstat=vector('numeric'),pvalue=vector('numeric'))
for(k in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[k])
if(uidlist[k]!=2305){
playmodel=glm(correct~incentive40+incentive70+incentive95,data=playframe,family=binomial())
playtest=conTest(playmodel,constraints=restrictmat,type="B",rhs=c(0,0,0))
addframe=data.frame(uid=uidlist[k],fstat=playtest$Ts,pvalue=playtest$pvalue)
}
#for subject 2305 who we check separately
if(uidlist[k]==2305){
addframe=data.frame(uid=uidlist[k],fstat=0.01,pvalue=0.99)
}
NIACindiframe=rbind(NIACindiframe,addframe)
}
sum(NIACindiframe$pvalue<0.05)/length(NIACindiframe$pvalue)

#NIAS and NIAC together
comboframe=data.frame(uid=vector('numeric'),NIASviol=vector('numeric'),NIACviol=vector('numeric'))
for(i in 1:length(uidlist)){
NIACplayframe=subset(NIACindiframe,uid==uidlist[i])
NIASplayframe=subset( NIASAframe,uid==uidlist[i])
addframe=data.frame(uid=uidlist[i],NIASviol=(sum((NIASplayframe$A_state_coefficient<0)*(NIASplayframe$Pvalue<0.05))>=1),NIACviol=(NIACplayframe$pvalue<0.05))
comboframe=rbind(comboframe,addframe)
}

comboframeclean=comboframe
NIASonly=100*sum((comboframeclean$NIASviol==1)*(comboframeclean$NIACviol==FALSE))/length(uidlist)
NIAConly=100*sum((comboframeclean$NIASviol==0)*(comboframeclean$NIACviol==TRUE))/length(uidlist)
both=100*sum((comboframeclean$NIASviol==1)*(comboframeclean$NIACviol==TRUE))/length(uidlist)
neither=100*sum((comboframeclean$NIASviol==0)*(comboframeclean$NIACviol==FALSE))/length(uidlist)
violatereport=data.frame(violate=c("NIAS only","NIAC only","Both","Neither"),percent=c(NIASonly,NIAConly,both,neither))

#************************************
stargazer(violatereport,summary=FALSE,rownames=FALSE)
#************************************